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SUMMARY 


An accurate and efficient method has been developed for the 
aerodynamic analysis of a series of arbitrary small geometry 
perturbations to a given baseline configuration. The method is 
appropriate for wings and wing- fuselage configurations in 
incompressible, potential flow. Mathematical formulations are 
presented for three computer programs that are employed. 

The first computer program is a conventional panel method. 
The surface of an arbitrary baseline configuration is represented 
by an array of panels on which piecewise constant source and quad- 
ratic doublet singularity densities are distributed. Based on 
the classical third identity of Green, the source strength cor- 
responding to prescribed normal velocity (Neumann) boundary 
conditions is calculated immediately from the local panel geome- 
try. The doublet strength is determined from the theoretically 
equivalent boundary condition that the internal perturbation 
potential vanish. A matrix of aerodynamic influence coefficients 
is calculated in order to establish a linear, algebraic system of 
boundary condition equations. The solution to the system of equa- 
tions provides the solution array of doublet strength at boundary 
condition control points. 

The second program calculates a matrix of partial deriva- 
tives. Each element of the matrix is the rate of change of 
doublet strength at one boundary condition control point with 
respect to the change in a surface coordinate of the paneled 
baseline geometry. The matrix is generated by a first-order 
mathematical expansion to the equations of the conventional panel 
method. The mathematical expansion is based upon analytical 
differentiation, not numerical. 

For each baseline configuration, the calculated quantities 
from the first two computer programs establish an input file for 
the third. The purpose of the third program is to calculate the 
surface pressure distribution and forces and moments for a series 
of perturbations to the baseline geometry. The third program is 
similar to the conventional panel method, with one major excep- 
tion: the solution doublet strengths are generated by linear 
extrapolation using the partial derivative matrix. This elimi- 
nates the two expensive steps of a panel method: calculating the 
aerodynamic influence coefficient matrix and solving a large 
system of linear, algebraic equations. The accuracy and effi- 
ciency of the perturbation analysis method are demonstrated by 
example calculations for a low aspect ratio, swept wing. 



INTRODUCTION 


In the past decade, several surface singularity panel 
methods have been developed for calculating the incompressible, 
inviscid pressure distribution on arbitrary three-dimensional 
configurations (references 1-9). A characteristic of these 
methods is good prediction accuracy for complex configuration 
geometries. Compared to an exact potential flow solution, the 
accuracy is limited only by the number of panels used to repre- 
sent the surface geometry. Nonetheless, surface panel methods 
are not widely used for the preliminary aerodynamic design of 
airplanes. The reason is computing cost. To analyze each change 
in geometry, it is necessary to compute a new, large matrix of 
aerodynamic influence coefficients and then solve a corresponding 
system of linear, algebraic equations for the solution singular- 
ity distribution. 

This report presents a method recently developed at 
McDonnell Aircraft Company (MCAIR) for the efficient analysis of 
a series of arbitrary small geometry perturbations to a baseline 
configuration. Designated the perturbation analysis method, it 
was developed under contract to NASA, Langley Research Center. 
The approach is to first calculate a conventional panel method 
solution for an arbitrary baseline configuration and then calcu- 
late a matrix consisting of the partial derivatives of 
singularity strength with respect to geometry coordinates. The 
baseline solution and derivative matrix are calculated one time 
only and then stored for repetitive use. For each geometry 
perturbation, the solution singularity distribution is 
constructed by multiplying the derivative matrix by a new 
right-hand side. This bypasses the two computationally expensive 
steps of a conventional panel method: calculating the influence 
coefficients and solving the large system of linear algebraic 
equations . 

The perturbation analysis method is especially appropriate 
for application to isolated wings and to wings on fuselages, 
because aerodynamic properties of wings tend to vary linearly 
with respect to thickness, camber, and twist. Unlike the exist- 
ing small disturbance methods, however, the baseline wing is not 
limited to a zero thickness plate. Furthermore, the perturbation 
analysis method can accurately predict the change in peak pres- 
sure coefficient with respect to changes in leading edge camber 
and radius of curvature. Even small changes in wing planform and 
fuselage geometry can be analyzed. 

To apply the perturbation analysis method, the following 
three computer programs are required (see figure 1): 

(a) MCAIR 3-D Subsonic Potential Flow Program (Version 2) 

(b) MCAIR 3-D Geometry Influence Coefficient Program 
(Version 1) 
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(c) MCAIR 3-D Perturbation Analysis Program (Version 1) 



MCAIR 3-D Subsonic 
Potential Flow Program 
(Version 2) 

MCAIR 3-D Geometry 
Influence Coefficient Program 
(Version 1) 

MCAIR 3-D Perturbation 
Analysis Program 
(Version 1) 

Input 

• Baseline Geometry 
(x, y, z)j 

• Baseline Geometry 
(x, y, z)j 

• Baseline Singularity 
Distribution n\ 

• Baseline Properties 
(x, y, z)j, n \ 

• Derivative Matrix 

(3 At j/9xj, 3^t j/3yj, 9pj/9zj) 

• Geometry Perturbation 
(Ax, Ay, Az)j 

Approach 

• Conventional Panel 
Method 

• First-Order Expansion 
to Panel Method 
Formulation 

• Linear Extrapolation 
A^j = 2 [(8/tj/3xj) Axj 

+ (9pj/9yj) Ayj + (9 jLt j/9zj ) Az ; ] 

Output 

• Baseline Singularity 
Distribution p j 

• Baseline Aerodynamic 
Properties 

• Derivative Matrix 

(9/V 9x j< dMi/^Vj. 9 At j/Dzj) 

• Aerodynamic Properties 
of Perturbed Geometry 


Figure 1. The Three Computer Programs for the Perturbation Analysis Method 


The first program (a) is an improved version of the conven- 
tional surface panel method of reference 9. The improvements, 
developed as part of the 1930 MCAIR Independent Research and 
Development Program, include a subpanel formulation that enforces 
doublet continuity at panel edges. The mathematical formulation 
for the program is presented in this report. For the 
conventional calculation of pressures, forces, and moments on an 
arbitrary configuration, only the first program is used. 

For the perturbation analysis method, the second program (b) 
employs an output solution file from the first program plus a 
differential mathematical formulation to calculate the partial 
derivative matrix. For each baseline configuration, the first 
two programs generate an input file for the third program (c). 
At any future date, the third program can be repeatedly executed 
at low computing cost in order to predict pressures, forces, and 
moments corresponding to arbitrary small geometry perturbations. 
Both the second and third programs were developed under .NASA 
Contract NAS1-16156, and the mathematical formulation is 
presented in this report. 
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In addition to documenting the formulation of the three com- 
puter programs, the purpose of this report is to present example 
solutions that demonstrate prediction accuracy and computational 
efficiency. A review of classical surface singularity theory for 
Laplace's equation and the application to panel methods is avail- 
able in reference 8 and is not repeated here. 

Use of trade names or names of manufacturers in this report 
does not constitute an official endorsement of such products or 
manufacturers, either expressed or implied, by the National 
Aeronautics and Space Administration. 
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DEFINITION OF THE ANALYSIS FLOW PROBLEM AND 
ANALYTICAL RESTRICTIONS 


Consider any number of bodies immersed in an incompressible, 
potential flow field (figure 2). A body must have finite thick- 
ness, and no point on the surface of one body can contact the 
surface of another. The fluid velocity is measured with respect 
to an arbitrarily selected Cartesian coordinate system (x, y, z) . 
It is assumed that either (1) the fluid is unbounded and would 
have uniform velocity if the bodies were removed or (2) the 

entire field is bounded by (contained within) one additional 
body. Suppose that the component of flow velocity normal to the 
surface is known at every body surface point in the field. In 
the analysis or Neumann problem, the objective is to solve for 
the two tangential velocity components. For steady state flow, 
Bernoulli's equation can be used to convert velocity to surface 
pressure, and subsequent pressure integration can provide body 
forces and moments . 




Ring Wing 


Figure 2. Bodies Immersed in a Flowfield 
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Sharp edges and points on a body are permissible and are 
designated “ridges" and "isolated points," respectively (figure 
3 ) . Formally, a ridge is a path on a body surface for which 
(1) the path is continuous, (2) the path direction is continuous, 
and (3) the body surface unit normal vector is generally discon- 
tinuous across the path. Examples of ridges are wing trailing 
edges and wing- fuselage intersections. An isolated point is a 
point where (1) every body surface path through the point has a 
discontinuity in direction and/or (2) three or more ridge end- 
points coincide. Examples of isolated points are the vertex of a 
cone and trailing edge point of a wing tip. Typically, a ridge 
will begin and end at isolated points. However, a ridge can be a 
closed loop with no isolated points, in which case the selection 
of an endpoint is arbitrary. 



Figure 3. Illustration of Isolated Points and Ridges 


In real flows, a thin viscous wake is convected from sharp 
edges by the fluid. In accordance with classical lifting surface 
theory, it is assumed that the effect of a wake on the flow field 
can be analytically simulated by introducing one or more trailing 
vortex sheets, herein designated "wake segments" (figure 4). It 
is also assumed that the geometry of each wake segment is known a 
priori. The four edges of a wake segment are designated the 
"forward edge," the "side edges," and the "rear edge." The side 
edges should be approximately aligned with the local flow 
direction. Restrictions on wake segments are presented below: 

(1) Every point on both surfaces of a wake segment must be 
wetted by the potential flow, 

(2) The length of each of the four edges of a wake segment 
must be finite (neither zero nor infinite). 
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(3) The classical Kutta condition is imposed on the forward 
edge, which must entirely coincide with a ridge of a 
body. If the ridge is longer than the forward edge, 
one or more additional wake segments must be attached 
to the ridge. In general, either no portion of a ridge 
or else the entire length must coincide with wake 
forward edges, 

(4) A forward edge cannot coincide with any other wake 
segment edges, 

(5) It is permissible for any portion of a side edge to be 
exposed to the potential flow, to coincide with another 
side edge (of even the same wake segment) , or to coin- 
cide with a body surface, 

(6) The rear edge, which corresponds to the starting vortex 
of classical lifting surface theory, must be fully 
exposed to the potential flow. The edge must not 
contact any part of any body or any wake segment. 


The computer programs described in this report are dedicated 
to solving the analysis problem defined above. The mathematical 
solution formulation is presented in the next section. 
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MATHEMATICAL FORMULATION 


The MCAIR 3-D Subsonic Potential Flow Program (Version 2) is 
a surface panel method that is based on the combined source- 
doublet distribution of the classical third identity of Green 
(reference 10). The advantages of this combined source-doublet 
distribution for a surface panel method are well documented 
(reference 8). The mathematical formulation employs a constant 
source distribution and a quadratic doublet distribution on each 
panel, where the solution singularity strengths are determined by 
satisfying indirect internal perturbation potential boundary con- 
ditions (reference 4). The flow velocity on each panel is then 
established from local velocity-singularity relationships associ- 
ated with Green's third identity, instead of the direct summation 
of the influences of the singularities on all the panels. 

The formulation to the MCAIR 3-D Geometry Influence 
Coefficient Program (Version 1) is similar, with one major 
exception. In each step of the solution process, the derivative 
of each quantity with respect to geometric coordinates is 
established, rather than the quantity itself. This process 
culminates in the generation of the matrix of partial derivatives 
of singularity strength with respect to geometry perturbations. 

The details of the mathematical formulation are presented in 
the following subsections. In each subsection for which mathe- 
matical formulas are presented, the first portion of the 
subsection presents the equations that are appropriate for the 
MCAIR 3-D Subsonic Potential Flow Program (Version 2), and the 
second portion presents the analogous differential equation for 
the MCAIR 3-D Geometry Influence Coefficient Program (Version 1). 
The third computer program [MCAIR 3-0 Perturbation Analysis 
Program (Version 1)] is nearly identical to the first. The only 
significant difference is that implementation of the matrix of 
partial derivatives eliminates the need to calculate aerodynamic 
influence coefficients and to solve a large system of linear 
algebraic boundary condition equations. 

Division of a Body Into Segments 

For computational purposes, the surface of every body is 
considered to be divided into a set of singly-connected surfaces 
designated "body segments" (figure 5). Each segment is bounded 
by four corners and four edges. A zero length edge is 
permissible. However, no zero length edge can touch a ridge 
(except at a ridge endpoint that is also an isolated point). To 
be consistent with surface singularity theory, every point on the 
wetted surface of a body should be part of a segment. 

With the exception of the following restrictions, the 
division of a body surface into segments is arbitrary: 

(1) Only edges (and corners) of a body segment can contact 
a ridge. 
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Figure 5. Division of a Body into Segments 

(2) Only corners and zero length edges of a body segment 
can contact an isolated point or a ridge endpoint, 

(3) Only edges (and corners) of a body segment can contact 
the side edge of a wake. Furthermore, such contact 
must occur along the full length of the body segment 
edge, 

(4) If a portion of a body segment edge coincides with a 
flow symmetry plane and the method of images is 
employed, then the full length of the edge must coin- 
cide with the plane. 

It is permissible for two edges of one body segment to coincide. 
For example, a single segment could represent both the upper and 
lower surfaces of a wing, in which case two opposite edges of the 
segment would coincide at the wing trailing edge. 

The total number of body and wake segments in the field is 
designated ng. If the method of images is employed, the x-z 
plane must be a flow symmetry plane and only the segments for 
which y> 0 are included in the total ng. If the x-y plane is a 
second flow symmetry plane, then only the segments for which both 
y>0 and z >0 are included. Each of the ng seqments is assigned 
an index ig (l£ig^ng) where the ordering is completely 
arbitrary. Body and wake segments are respectively identified by 
the nomenclature STYPE(ig) = 1 and 2. 
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Each of the four corners and edges of a segment are identi- 
fied by the indices "icRN" and " ^EDg" (figure 6). For a body 
segment, the ordering is counterclockwise with respect to an 
observer in the flow field looking at the surface. For both body 
and wake segments, the local surface normal vector is positive if 
it points out of the plane of the page in figure 6. For a wake 
segment, ijsDG — 1 is the forward edge, 2 and 4 are the side 
edges, and 3 is the rear edge. If the forward edge of two adja- 
cent wake segments are attached to the same ridge, the adjacent 
side edges must be i^DG = 2 of one wake segment and igpQ = 4 of 
the other. 



'CRN “ 1 

Figure 6. Ordering of Corners and Edges on a Segment 


Each corner icRN ° f each segment ig is categorized by 
CTYPE ( icRN' is ) where 

1 if the corner does not coincide with an isolated 

CTYPE = point, 

2 if it does. 

Each edge i^DG °f each segment ig is categorized by both 
ETYPA ( ijEDG • f S ) and ETYPB ( iEDG • is ) where 
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ETYPA = 


1 if the length of the edge is finite, 

2 if the edge is a zero length body segment edge 

not on a ridge, and 

0 if the edge does not lie on a symmetry plane, 

ETYPB = 1 if the edge lies on the y = 0 symmetry plane, 

2 if the edge lies on the z = 0 symmetry plane. 

Illustrated in figure 7 are body segments on opposite sides 
of a ridge. One side is selected to be designated igjDE = 1 and 
the other isiDE = 2 * Although it is not consistent with analyti- 
cal theory, it is sometimes reasonable to ignore the body on one 
side of a ridge. For example, the surface of a wing tip is 
usually not represented in practice. However, if the forward 
edge of a wake is attached to a ridge, then both sides of the 
ridge must be represented by body segments. 


Ridge 



Figure 7. Body Segment Edges on a Ridge 
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The total number of ridges that lie along edges of the ng 
body and wake segments is defined as n R , where n R is zero if 
there are no ridges. Each ridge is assigned an index i R 

( 1 _< i R _< n R ) . Usually, the endpoint of a ridge is an isolated 

point. However, if the ridge is a closed contour on which there 
is no isolated point, then one point on the contour must be 
selected as the endpoint. Furthermore, if a ridge passes through 
a symmetry plane and the method of images is applied, then the 
point on the symmetry plane is designated a ridge endpoint. 

Division of a Segment Into Panels 

The geometry of each segment ig is represented by a grid of 
points, where the points are identified by rectangular indexing 
( i , j ) 7 see figure 8. The value of the limits m and n on the 
indices i and j are selected by the user and can vary from one 
segment to another. It is usually recommended that each point 

which lies on the edge of one segment should coincide with a 

point on the edge of an adjacent segment. However, discontinui- 
ties in point density and even geometric gaps are permitted 
between adjacent segment edges. 


(1 ,n - 1) 







(m-1 , n-1 ) 













z^Panel ^ 












(1.2) 








(1, D 

(2. 1) 






(m-1 , 1) 


I Corner Point index (i, j) 

* Panel Index (ip, jp) 

Figure 8. Schematic of a Segment Divided into Panels 

For consistency with the segment corner nomenclature of 
figure 6, it is stipulated that the values i CRN = {1; 2; 3; 4 } 

correspond to the points (i,j) = {(1,1); (m,l); (m,n); (l,n) }. 
Consequently, i-EDG = {li 2? 3; 4} corresponds to {j=l; i=m; j=n; 
i=l ). 
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The four points {(i,j); (i+l,j); (i+l,j+l); (i,j+l)} are the 
four corners of a geometric panel. Each panel is identified by 
the pair of indices ip and j p . The edges of a panel are formed 
by the four straight line segments that connect adjacent panel 
corners, where the line segments are not necessarily coplanar. 

A cumulative index k identifies individual points for the 
collective system of ng body and wake segments. The conversion 
from the point indexing (i,j) of segment ig to the cumulative 
index k is 


where 


W 


k = k Q (ig) + ( i-1) *n + j 


( la) 


(ig-D 


j S =1 


^ = 1 


[m( jg) -n ( jg) ] 


2 1 i s 1 (n s +l) 


(lb) 


The geometry of the entire configuration is defined by the follow- 
ing array of coordinates to be specified by the user: 


(x.y/zjfc where 1 _< k k 0 (ng+l) 


A cumulative panel index kp is related to the panel indexing 
(ip'jp) °f segment ig as follows: 


k P = k P {i S } + (ip-D • (n-1) + j p (2a) 
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where 


0 . 


•V 1 


(2b) 


k p (ig) 
o 


(ig-D 

l [m( jg) -1] • [n ( j g) -1] 
j s =l 


2 £ ig <_ (n g +l) 


Panel Properties 

Presented in this subsection are the formulas for calcula- 
ting the properties of a panel as a function of the coordinates 
of the four corner points. For each panel, a quadratic doublet 
distribution is established by a least squares fit to the doublet 
strength at several spline control points on the panel. Because 
the doublet strengths are initially unknown, a sensitivity matrix 
is generated for which each matrix element is the contribution of 
a unit strength doublet density at one spline control point to 
one of the six coefficients of the quadratic distribution. Also 
presented in this section is the procedure for generating a 
second sensitivity matrix, one that relates the change in quad- 
ratic doublet coefficients to panel geometry perturbations. 

Consider the panel of figure 9, for which the four nonco- 
planar corners k = {1,2, 3, 4 } respectively correspond to the 
points {( i , j ) , (i+l,j), (i+l,j+l), (i,j+l)} of figure 8. Define 
the array X-^ as the coordinate displacement between points 1 and 
k: 


'V 


’ x k " x l' 

* x 2k 

v — 

III 



y k - y l > 

. X 3k 


. Z k - Z l. 


where (x,y,z) refers to the global coordinate system of figure 2. 
The panel geometry is uniquely determined by nine quantities that 
are assumed to be known: (1 < i < 3, 2£k_^4). In accordance 
with equation (3), X^ k is identically zero for k = 1. The points 
designated k = {5, 6, 7, 8} in figure 9 are, respectively, the mid- 
points of the line segments connecting points {1 and 2, 2 and 3, 
3 and 4, 4 and 1}; panel center point k = 9 is the midpoint of 
the line segment between points 5 and 7, which is coincident with 
the midpoint of the line segment between points 6 and 8. It is 
easily proved that the five points k = {5, 6, 7, 8, 9} are coplanar. 
x ik for these five points is calculated from the following 
formula. 
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2 [X i(k-4) 


+ X . /, -y \ 
i (k-3) 


5 < k < 7 


X ik 


= < 


2 [X i(k-4) 


+ X il ] 


k = 8 


|[ x i 5 + X i? ] 


k = 9 



Corners (noncoplanar) 


Edge midpoints 


Center 


Figure 9. Definition of Point Index “k” for a Panel 


(4) 


The coordinates of the panel center are designated (x,y,z) 0 . 

(x,y,z) o = (x^,y^/Z^) + ^ X 19 'X 29 ' X 39 ) (5) 

The panel area is calculated by the magnitude of a vector cross- 
product . 

AREA = | <X 17 - X i5 > X (X i8 - X i6 ) | ( 6) 
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The nine points k are designated "spline control points." 
Suppose that an arbitrary value of doublet strength (y^) is 
assigned to each spline control point. It is desired to estab- 
lish a quadratic distribution that will exactly match y-^ at the 
panel center and edge midpoints and approximately match y ^ at the 
panel corners. This is accomplished by applying the method of 
Appendix A, which takes advantage of the fact that points k = 5, 
7, and 9 are colinear and points k = 6, 8, and 9 are colinear. 
First, it is necessary to introduce a new set of indices and a 
local panel coordinate system that are consistent with Appendix 
A. 


As presented in figure 10, the index £ is introduced to 
designate the points 1< k £ 9 in a different order. The nomencla- 
ture £(k) and k(£) respectively denote the conversion of k^.to 
£ and £ to k. Define the array (equivalent to a vector P^) 

as the coordinate displacement of each point £ from the panel 
center: 


P i£ " P £ “ [X i,k (£) X i,k (1) ] 


(1,1) < ( i , £ ) < (3,9) (7) 
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A local ( ti, £) Cartesian coordinate system with origin at 
the panel^ center is constructed (figure 11) such that the unit 
vectors (e^e^e^) are 







( 8 ) 


e = (P 2 x P 3 ) / 1 (P 2 x P 3 ) 


( 9 ) 


e 

n 



x e ? ) 


( 10 ) 
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The vector e ^ is designated the panel unit normal vector. Equa- 
tions (8) - (10) are equivalent to a rotation matrix a^j, where 




-> 

e. 


e 

5 


X 

e 

- = [a. . ] H 

-* 

e 

9 

l j 

y 




e 


e 

^ > 


z 


( 11 ) 


The angle 0 of figure 11 is 


e = \ cos - 1 [(? 2 -P 3 )/(|? 2 | • \ $ 3 \ ) ] 


(12) 


Define the distance between the panel center and point JL as 


0 


£ = 1 


s 


£ ~ 



£ = 2,3 - 
£ = 4,5 . 


(13) 


The coordinates of the panel corners in the local (5 # ri#C) 
coordinate system are 


J 


? £ = P £ 


-> 9 

e „ 


->- -> 
p • e 
£ n 


6 < £ < 9 


(14) 


A coplanar adjustment is performed by setting 


= 0 (6 £_ £ <_ 9 ) 


(15) 
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With the aid of the above defined quantities, the method of 
Appendix A will establish the quadratic doublet distribution in 
the following form: 


yU,n) = ^(i) + + e 2 n + + ^ 4 ^ 


+ 3 5 * (n 2 - 5 2 tan 2 0) 


(16) 


where 


B i [b i2, ‘ y kU) 1 + <S i5* £ ^ 6 [b 5fi,‘ p k(2,) ] 


(1 < i < 5) (17) 


Each quadratic coefficient is a linear function of y^. 
The sensitivity matrix elements bj_£ [(1,1) £ (i,&) £(5,5)] and b 5 £ 
(6 ££.£ 9) are calculated independently of y^; in fact, the calcu- 
lations depend only on the nine quantities Xq] c ,X 2 fc,X 3 ] C (2£ k£ 4). 
For k = 1, the quantity Xjj^ is always zero and therefore does not 
affect the calculations. 

The preceding discussion addressed the procedure for estab- 
lishing the quadratic doublet distribution corresponding to known 
values Xj_fc and arbitrary unknown values y^. Now suppose that 
each is assigned a fixed, known value. The method for deter- 
mining the first-order change in quadratic doublet distribution 
induced by geometry perturbations (dX^^) is presented below. 

Define the array "v" as follows: 


^ V l' v 2 ' V 3' V 4 ' V 5 ' V 6 ' V 7' V 8' v 9^ 


- {X 12 ,X 22 ,X 32 ,X 13' X 23' X 33' X 14' X 24 ' X 34 } 


(18) 
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I 


The purpose behind the introduction of single-subscripted array 
"v" is to simplify the indexing. Then dX-j^ for l£k£9 can be 
expressed as 


3 

dX.. = E DXI. .*dv 

ik j = 1 kj [i+3(d-1)] 


1 < i < 3 


(19) 


where 


DXI. . s 
k3 


r 




< 


k ( j+1) 


2 [DXI (k-4) j + DXI (k-3)j ] 

i [DXI (k-4)j + DXI lj ] 
|[DXI 5j + OXI 7j ] 


if k = 1 


if k = (2,3,4) 


if k = (5,6,7) > (20) 


if k = 8 


if k = 9 


J 


Equations (19) and (20) give the rate of change of the 27 values 
Xfk with respect to each of the nine members of the v-array. 

The derivative of Pj_ ^ can be expressed as 


dP i£ = ^ DPI Ij‘ dV [i+3(j-l)l (1,1) i (i.t) i (3,9) (21) 


where 


DPI *j 


[DXI 


k («,) , j 


- DXI 


k ( 1) , j 


] 


( 22 ) 
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By applying the calculated values DXI^j and DPI^j to the 
derivatives of equations (8) - (14), the following quantities are 
calculated for 1 < m < 9: 


3 a . . 
ii 


30 


3 s . 


3 v ' 3v ' 3 v 
m mm 


l for (1,1,2) < (i , j , & ) < (3,3,5) 


and 


35, 


3n 


3 v ' 3 v 
m m 


for 6 < Z < 9 


The method of Appendix A will provide calculated values for 
93i/90» 3|3i/3s^ for 2 ££,<_ 5, and (30^/31^, Sj^q/Sq^) for 6 9. 
The change in quadratic doublet coefficient (d|3q) generated by 
panel geometry perturbations is 



All of the partial derivatives of equation (23) are calculated 
quantities. Consequently, the multiplications and summations of 
equation (23) can be performed to provide ^^i/^ v m for 
(1,1) < (i,m) < (5,9) . 

Now consider arbitrary simultaneous perturbations to both 
panel geometry (dv m ) and doublet strength at the spline control 
points (dy^) . Then the change in quadratic doublet coefficient 
(d(3q) can be predicted by the following formula. 
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( 24 ) 


de 


9 

I 

m=l 


17 3 3 . \ 

' 

L 1 

dv 

\ 3v 

m 

L\ m/ 

. 




+ S i5‘,,' 6 tb 5t' d “k(t) 


All the coefficients of equation (24) are known quantities, calcu- 
lated by the procedure described above. In the same sense that 
the linear relationship between and of equation (17) is 
appropriate for a conventional panel method, the linear relation- 
ship between dg-j_ and dv m of equation (24) is appropriate for the 
analysis of small geometry perturbations. 

The formulation presented above in this section adjusts each 
noncoplanar panel in order to form a planar quadrilateral. This 
will typically generate small gaps between the edges of adjacent 
panels. Furthermore, the quadratic doublet distributions will 
typically be discontinuous across the edge. These small geo- 
metric and singularity discontinuities are usually insignificant, 
except for panels that are adjacent to a ridge. An alternate 
formulation that eliminates these discontinuities at the expense 
of increased complexity is presented in Appendix B. There each 
panel is divided into eight triangular subpanels. The alternate 
formulation is always used for each panel that is adjacent to a 
ridge. For the remaining panels, the formulation of this section 
is normally applied. 

Boundary Condition Control Points 

A set of control points on body and wake segments is speci- 
fied as locations where flow boundary conditions are to be 
satisfied. The center of every panel of every body segment is 
always a control point. In addition, certain body segment points 
along ridges and certain points on the forward edge of each wake 
segment are boundary condition control points. 

Every control point is assigned an index k^p. The range of 
indices for body or wake segment ig is designated 


"CP 


(ig) 


+ 1 


CP - 


: CP (l s +1) 

O 


where k C p Q (i s ) is the total number of control points on the 
segments preceding ig. 
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The boundary condition control points on the panels of body 
segment ig are established in the following orders 

(1) The center of each panel (ip,jp) is a control point, 
and the index is 


k CP k CP + (i p -1) * (n-1) + jp 


(2) The midpoint of each panel edge that adjoins a ridge is 
a control point. The control point ordering proceeds 
counterclockwise around the segment from iEDG = 1 
through i ED G = 4 * Exception: There is an option to be 
discussed later for which all the panel edge control 
points are eliminated on the side of the ridge desig- 
nated by the user as isiDE = 2. 

(3) If a corner of a body segment coincides with an 

isolated point, there is normally a control point at 
the corner, and the control point indexing increases 

with respect to increasing corner index icRN • The 

control point will exist if CTYPE( i^p^, ig) = 2 and 
ETYPA( i C RN~l/ ig) = 1 [where (icRN -1 ^ is to k e replaced 
by 4 iff icRN is 1] • 

Each boundary condition control point of a body segment is 
located near but not on the local panel (or subpanel) surface. 

The locatio n is on the interior side of the panel at a distance 
/panel area • 10 _ ® from the surface. The control points on panel 
edges and corners, (2) and (3) above, are retracted slightly in 

the general direction of the panel center. The approximate 

distance retracted is a fraction [e(ig)] of the distance to the 
panel center, where £(ig) is a value specified by the user 

[recommended value for body segments: e (ig) = 0.02]. 

Along the forward edge of each wake segment ig, control 

points are established in pairs. For each pair, one point (k c p) 
is slightly below the local wake panel (or subpanel) and the 
second point (kcp+1) is slightly above. ^ The line integral of 

velocity between the two control points (/V* ds) will be set equal 
to zero, thereby satisfying one Kutta condition equation for each 
pair of control points. For wake segment i s , each pair of boun- 
dary condition control points is established in the following 
order: 


24 



( 1 ) 


There is one pair of control points centered at each 
panel edge midpoint along the segment forward edge 
( iEDG = 1) • The two control point indices are 


k cp = k cp (i s> + + n 2 > 


(2) There is one pair of control points centered at segment 
corner icRN = 1 if CTYPE(l,ig) = 2. Furthermore, there 
is a pair at Icrn = 2 if CTYPE(2,i g ) = 2. 

Each pair of wake control points is retracted from the panel 
edge or corner in the same manner described above for body 
control points. For wake segments, however, the recommended 

value for the retraction factor e(ig) is generally not the value 
of 0.02 that is appropriate for body control points. Instead, 
the guidelines in the following paragraph should be employed. 

Figure 12 is a cross-section view in which the forward edge 
of wake segment i g is attached to a ridge of a body. Typically, 

the ridge will be the trailing edge of a wing. The lengths As^, 

As B , and As c are the dimensions in the cross-section view of the 
three local panels that adjoin the ridge. The distance between 
the two wake control points (k^-p and k^-p+l) is approximately 

equal to the retraction distance from the ridge. It is recom- 
mended that the value e(ig) be specified such that 


eUg) 


0.02 x 


/ AS A + tS B 
\ AS C 


J'or each pair of wake control points, the unit direction vector 
Wn is specified by the user. In typical practice, Wjvj is normal 
to the local bisector of a wing trailing edge. 

For each body control point and each wake control point, the 
calculated coordinates are designated (x,y,z)^ cp . Corresponding 
to the discussion above, these coordinates are a function of the 
locations of the four corners of the local panel and, in the case 
of a wake control point, a function of the local vector If 
the subscripts 1, 2, 3 and 4 are used to denote the four panel 
corners, then the first order change in control point location 
induced by perturbations to both the panel geometry and W N can be 
expressed as 
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Figure 12. Pair of Wake Control Points k cp and k cp + 1 


d (x,y , z) , = + a-d(x,y,z) +b-d(x,y,z) 2 

k CP 

+ c-d(x,y,z) 3 +d-d(x,y,z) 4 (25) 

+ e • d (W ,W ,W N ) 
x y z 

where the coefficients (a,b,c,d, e) are calculated quantities. 
For example, if the control point is at a panel center, then the 
coefficients a,b,c and d are each 0.25. The value of e is always 
zero if the control point is on a body segment. For each pair of 
wake control points, the effect of distance between the points is 
considered to be computationally insignificant, and the change in 
distance is therefore ignored in the perturbation analysis 

method. 

In addition to the boundary condition control point indexing 
kcp described above, a compressed indexing designated k' is 
defined. The ordering of k' is the same as kcp, except that for 

each pair of wake control points k^p and kcp+1, there is only one 

value k^ p . Therefore, the total number of indices k^p is equal 
to the number of body boundary condition control points plus one- 
half the number of wake points. For each k^p there corresponds a 
local unknown doublet strength that is to be determined by 
satisfying a local flow boundary condition. 
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Doublet Interpolation Between Neighboring Control Points 


The section entitled "Panel Properties" described the proce- 
dure for establishing a quadratic doublet distribution for a 
panel on the basis of nine spline control points. For each panel 
on a body segment, one or more of the nine points is a boundary 
condition control point, and the local doublet strength is an 
independent unknown that is to be determined by solving a large 
system of linear, algebraic boundary condition equations. For 
the other points, however, the strength is established by inter- 
polation between neighboring boundary condition control points. 
Interpolation is also performed along the forward edge of each 
wake segment, where there is an independent unknown doublet 
strength for each pair of boundary condition control points. The 
interpolation procedure for body and wake segments is presented 
in this section. 


Figure 13 is a schematic of all the doublet spline control 
points of a segment. Each point is identified by the indices 
(iy,3y). For a wa ^ e segment, the doublet strength at each point 
( i , j y ) is specified to be identical to the strength at the 
forward edge point (iy , j y =l). This corresponds to the conven- 
tional wake condition that the streamwise doublet gradient is 
zero. For each wake segment, the variable index jy is therefore 
replaced by the constant jy = 1. Cumulative indexing is defined 
for the doublet spline control points of all segments 


(i -1) • ( 2n-l) + j if STYPE ( i c ) = 1 

y y fa 


k 

p 



4 - 


where ky Q (ig) is 
preceding ig. If 
and the method of 
k y is assigned the 


r (26) 


i if STYPE (i s ) = 2 

> 

the total number of ky's on the segments 
there is one or more flow planes of symmetry 
images is used, then the mirror image of point 
same index ky. 


The doublet strength for each ky is to be a weighted average 
of the doublet strengths at np neighboring ky'S, where each neigh- 
boring ky must correspond to a boundary condition control point 
(or to a pair of wake boundary condition control points). The 

neighboring locations are identified by the nomenclature ky^, 
where i = 1,2,3, ...,np. The objective is to determine the points 
k y ^ and weightings c_-j_ such that 


n T 


pOO 


= £ c . * p (k ) 


i=l 


p 


(27) 
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(2n — 1 ) 



(2m - 1) 


Figure 13. Doublet Spline Control Points on a Segment 


If ky corresponds to a boundary condition control point, the 
weighting is trivial? i.e., np = 1, ky^ = ky, and c^ = 1. In 

general, quadradic interpolation is required to establish the 
proper weighting. Surface interpolation is applied to the points 
on body segments. Then for body and wake points along ridges, 
path interpolation is applied. 

Figures 14 and 15 are surface interpolation schematics of a 
point ky on a body segment and the neighboring points i 
(l£i^np). Figure 14 applies if ky is a panel corner, i.e., if 
the values iy and j y are both odd. Figure 15 applies if iy is 
even and jy is odd. Interchanging iy and jy in figure 15 corre- 
sponds to odd iy and even jy . If the separation between ky and 
every edge of the segment on which ky lies is greater than one 
panel, then all of the np neighboring control points are at panel 
centers. Otherwise, some of the np points will either be on the 
opposite side of the segment edge or else not be panel centers. 
The former applies if the edge is neither part of a ridge nor a 
zero length edge at an isolated point, and it is the responsi- 
bility of the user to select the values on the opposite side of 
the edge. Otherwise, the indexing is handled automatically, and 
a single edge or corner control point can be identified by more 
than one index i. 

It is noteworthy that there are unusual cases for which the 
eight point interpolation of Figure 15 is less accurate than 
twelve point interpolation. Therefore, the consistent use of 
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Figure 14. Local Indexing at a Panel Corner Point 
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Figure 15. Local Indexing at a Panel Edge Point 




twelve point interpolation in the formulation is being 
considered. 

The quadratic surface interpolation is not applied to the 
conventional doublet strength y but, rather, to a modified 
doublet quantity designated u. As illustrated in figure 16, this 
accommodates the situation in which a body segment edge that is 
not a ridge coincides with a side edge of a wake. In the 
vicinity of point ky , u is defined to be identical to y . How- 
ever, on the opposite side of the segment edge, u is adjusted by 
+ Ay, where Ay is the local wake doublet strength. Therefore, the 
quantity u corresponds to a continuous function with continuous 
gradient on the body surface. This makes u a suitable quantity 
for quadratic interpolation, whereas interpolation on y would be 
inappropriate across the edge. 




Figure 16. Doublet Discontinuity at Wake-Fuselage Intersection 


Before the surface interpolation is initiated, it is neces- 
sary to establish a local (£/H/ l ) -coordinate system such that 
the origin is at ky and the £-axis is locally normal to the body 
surface. The vector Pj_ is defined for l£i<np. 
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( 28 ) 


where the positive sign for y and z is selected unless ky^ is on 
the opposite side of the y = 0 and z = 0 symmetry planes, respec- 
tively. Vector Pi is the displacement between points ky^ and ky. 
The unit vector e ^ is a weighted average of local panel normal 
vectors. 


n. 


n , + n , + n 
x' — ^ ' 


'? _ 


kp(i) 


1=1 / 


2 2 2 
P Z +P. Z +P. 
ll i2 i3 


(29) 


where the constant c accounts for the fact that the right-hand 
side is not normalized, where kp(i) is the panel on which point 

k yi lies, where (n x ,ny,n z ) is the unit normal vector of panel 
kp(i) in the ( x, y, z) -coordinate system, and where the summation 
is to include only values i that are panel centers less than one 
panel distant from ky . Once again, the sign selection depends 
upon whether point i is on the opposite side of a symmetry^ plane . 
e^ is defined^ by selecting any^ unit ^vector normal to e £. It 
follows that e n is defined as er x e£. Each point i is then 
assigned coordinates in the ( £ , n , C ) -system. 



(30) 
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In order to approximately account for surface curvatuve, S j_ and 
n i are both magnified by the factor 


L(d + n-t + cf)/(d + 


)] 


1/2 


Function u(£ ,n ) is modeled as a quadratic 


u(5 ,n) = 3 -]_ + + 


e 3 n 


+ 


+ 


V 


(31) 


Designating ui as the value of u at point ky^ (1£ ifno), 
quadratic interpolation is performed by minimizing function E. 


E = 


n D . 
Z W. ‘ 

i-1 1 


(3 


1 + S 2 5 i + B 3 n i + S 4 E i + e 5 E i n i + s 6 ri i 2 * V 


(32) 


where Wj_ 2 is 1.00 unless ky^ is separated from ky by less than 
one panel, in which case W^ 2 = 10^. E is minimized with respect 
to the unknowns The resulting system of linear 
equations is manipulated by standard matrix inversion and matrix 
multiplication techniques to provide the array c-^ where 


6 


1 


n T 


i=l 


c i u i 


(33) 
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Since the origin of the ( 4 , rp 5 ) -coordinate system is at point ky , 
equation (31) implies that 


y (k ) 
y 


(34) 


Combining equations (33) and (34), 


y<V 


E 

i=l 


c . u . 


1 1 


(35) 


Equation (35) can be transformed into the desired form of 
equation (27) by expressing u^ as y(ky ^ ) + Ayq + Ay 2 + •••/ where 

each Ayij corresponds to the side edge of a wake segment intersect- 
ing the body surface between points ky and ky^. The value Ayj is 
equal to the doublet strength at a boundary condition control 
point on a corner of the wake segment. Therefore, equation (35) 
becomes 


y (k ) 
y 


'W 

1 

i=l 


c . 
1 


(k ) 


(36) 


where i >np designates the wake segment corners associated with 
Ayj [specifically Ayj -u(ky + .)]. If the limit np is increased 
to include rtyj, then equation (36) is identical to equation (27). 

The quadratic surface interpolation scheme described above 
is not applied to the points ky that lie on a ridge. Instead, 
the interpolation is based upon distance along the ridge. It is 
the responsibility of the user to identify the ordering of the 
segment index and segment edge index (is and isDG ) along each 
ridge. Then the pathwise interpolation method described below is 
automatically employed for the points on the ridge. 

Illustrated in figure 17 are four panels along a ridge. The 
panels can be either body panels or else be wake panels along a 
wake forward edge. The four panels need not belong to the same 
segment. The objective is to establish weightings cj_ for 
equation (27) where point k y is at the panel corner depicted in 
figure 17, index i refers to the four panel edge midpoints, and 
np = 4. If the point ky is one panel width from a ridge 
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endpoint, then one of the four midpoints i either is replaced by 
the ridge endpoint or else is the first panel edge midpoint past 
the ridge endpoint. This respectively depends upon whether the 
endpoint is or is not an isolated point. Similarly, if ky is at 
the ridge endpoint, then either the interpolation is trivial 
because ky is also a boundary condition control point or else two 
of the points i are past the endpoint. 



The interpolation is based on the distance s in figure 17. 
A parabola is established that matches the doublet strength 
exactly at points i = (2,3) and approximately at points i = 

(1,4). Define s' such that the origin is midway between i = 2 
and 3 . 


s ' = s 



+ s 3 ) 


(37) 


Exact doublet matching at i = 2 and 3 implies that the quadratic 
11 ( 3 ' ) is of the form 


V (s') 





+ 


y 3 ) 


(38) 
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where a is a constant. A weighted square error in matching 
doublet strength at i = 1 and 4 is defined as 


E = (S-^ 1 ) 4 * [u ( s ^ ' ) - V^] 2 + (s 4 ' ) 4 - [\i (s 4 ' ) - u 4 ] 2 (39) 

The constant a is evaluated by substituting equation (38) into 
(39) and minimizing E with respect to a. Then the value s' 
corresponding to point Ky is applied to equation (38), which 
reduces to the desired form of equation (27). 

Heretofore, it has been assumed that the points i = 1,2, 3,4 
in figure 17 are boundary condition control points. However, for 
each ridge there is an option which eliminates the need for 
boundary condition control points on the side designated by the 
user to be "i-siDE = 2" (see figure 7). For each point ky on 
-*-SIDE = 2 (i.e., for each panel corner and each panel edge mid- 

point), the doublet strength is established directly from the 
theoretical condition that the net local doublet strength from 
-i-SIDE = 1 plus isiDE = 2 plus wake (if present) must be zero. 
This option is usually employed because it saves computing 
effort . 

For every spline control point ky , the above formulation of 
this section determines a set of points ky ^ and weightings ci 
such that ij(ky) is established by equation (27). Each ky^ is 
equivalent to one boundary condition control point index k^p. 
Using the nomenclature k^p Hk ( i,p(ky^) to denote the transforma- 
tion from index ky^ to k^p 1 , equation (27) is cast in the final 
form 


P(k u ) 


E 

i=l 


c i* ,i(k CP. ) 

l 


(40) 


It is noteworthy that each weighting c^ is a function of the 
panelled geometry coordinates but not the doublet strengths. In 
the perturbation analysis method, the first-order change of cj_ 
with respect to geometry perturbations could be included. How- 
ever, the change is considered to be insignificant for typical 
applications of the perturbation analysis method and is neglected 
in order to reduce both formulation complexity and computing 
effort. Therefore, the differential form of equation (40) for 
the perturbation analysis method is 
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(41) 


dy (k ) 


E C . -dy (k' p ) 
i=l 1 ^ i 


Singularity Distributions on a Panel 

The effect of each panel on the flow is represented by a con- 
stant source distribution and a quadratic doublet distribution. 
Because the panel method of this report is based upon the inte- 
gral surface singularity theory associated with Green's third 
identity (reference 10), the source strength can be determined 
directly from the local panel geometry, as will be described 
below. In order to determine the doublet distribution, it is 
necessary to establish and solve a system of boundary condition 
equations, which is the subject of a following subsection. Here 
it will be appropriate to establish the functional relationship 
between the unknown doublet strength at each control point k' 
and the quadratic doublet distribution on each panel. Ct 

For each body panel, the constant source strength correspond- 
ing to Green's third identity is 


CT 



(42) 


where V^ p is the prescribed normal velocity boundary condition at 
the panel center. If the panel is divided into subpanels 
(Appendix B), then e^ refers to the subpanel normal vector. For 
wake panels, the source strength o is identically zero. 

The quadratic doublet distribution on each panel is related 
to the doublet strength at the local spline control points by 
equations (16) and (17) [analogous equations (B7) and (Bll) apply 
for subpanels]. Each of the nine spline control points in equa- 
tion (17) corresponds to one value k^ . Consequently, equation 
(40) can be combined with equation (17) to establish the relation- 
ship between each quadratic coefficient Si of equation (16) and 
the unknown doublet strength for 
ience, a more conventional set 
( 1 f. i 1 6 ) is defined as follows: 


MS,n) = a i + a 2 ^ + a 3 n 


each point k^p. For conven- 
o£ quadratic coefficients 


+ a 4 5 2 + a 5 £n + a g n 2 (43) 
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For each panel, equations (16), (17), (40), and (43) are manipu- 

lated to establish quantities MAX, k^ p , and d^j where 


a . = 


MAX 
Z d 

j = l 


ID 


y ^ k Cp 


1 < i < 6 


(44) 


The calculated quantity d^j depends upon the defining geometry 
coordinates but not upon the doublet distribution. In accordance 
with equation (44), each quadratic coefficient varies linearly 
with respect to each unknown doublet strength V(^£ P )» If the 
panel is divided into subpanels, there is one set of equations 
(43) and (44) per subpanel. 

For the perturbation analysis method, it is necessary to 
establish the first-order relationship between geometry perturba- 
tions and the resulting perturbations to source strength 0 and 
quadratic coefficients aj_. The source perturbation is estab- 
lished from equation (42), 


da 



(45) 


where the relationship between de^ and perturbations to the panel 
coordinates was established in the subsection entitled "Panel 
Properties" (or Appendix B for subpanels). For the doublet 
distribution, equations (24) and (41) (or analogous equations 
(B20) and (41) for subpanels) establish the first-order expres- 
sion for perturbations to the quadratic coefficients correspond- 
ing to perturbations to both panel geometry and to the unknown 


doublet strength at each point k^ p . Recall that the geometry of 

each panel is defined by the array v-; [equation (18)]. Then for 

J ^ 9 a i 

each panel (or subpanel), the quantities ( K Q ■■ ) and (~r — ^) can be 

O V j C7 V j 

calculated, where 


da 


9 


Z 

j=l 


3a . , 

) • dv . 

3v. D 


(46) 
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/ i \ 

da . = El I • dv . + 

1 j=i\ 3v j/ ^ 


and where dij is the calculated coefficient associated with 
equation (44). 

Potential Influence Functions 

Consider the planar quadrilateral panel of figure 11. With 
respect to the local ( £, n, 0 -coordinate system, the panel geo- 
metry is defined by the coordinates of the panel corners [(£]£,%) 
l^k^4]. Now consider the seven following singularity distribu- 
tions on the panel: 


(1) 

y (£# n) 

= 

1 

Constant Doublet 

(2) 

y ( C/ n) 

= 

K 

First-Order Doublet 

(3) 

y (S / n ) 

= 

n 

First-Order Doublet 

(4) 

y (S ' n) 

= 

E 2 

Second-Order Doublet 

(5) 

y (Z, n) 

= 


Second-Order Doublet 

(6) 

y^»n) 

= 

n 2 

Second-Order Doublet 

(7) 

0 • n ) 

= 

1 

Constant Source 


Define <p i as the potential at arbitrary field point (£,n,C) 
induced by the ith singularity distribution above. Mathematical 
formulas are available for calculating <j>i as a function of eleven 
quantities: n l» ^2» ^2* ^3» ^3# ^4» 14,? , r) , £} (see reference 2). 
The formulas are also appropriate for a triangular subpanel 
(figure B-2) if the subpanel is considered to be a quadrilateral 
for which the length of one side is zero. 

Analytical differentiation of the formulas for 4>i estab- 
lishes formulas for each of the following partial derivatives: 


MAX 

I 

j=l 


d ij' dl,(k ip. 


(47) 


3 ^' IT' 



(1,1) < (i,k) < (7,4) 
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For application to the perturbation analysis method, the first- 
order change in <J>i induced by perturbations to both panel 
geometry and field point location is expressed as 



Boundary Condition Satisfaction 


For each control point k^p, there is one boundary condition 
equation to be satisfied. On bodies, the prescribed normal velo- 
city condition is represented by an equivalent internal condition 
that the net potential induced by all the singularities is zero. 
On wakes, the Kutta condition is satisfied by requiring that the 
average flow velocity component vanish between each pair of wake 
control points k^p and k^p+l; this is equivalent to specifying 
that the difference in total potential between the two points is 
to be equal to the local wake doublet strength. 

Each boundary condition equation is established by summing 
the contributions from each panel (or subpanel). Using equation 
(43) and the definition of the influence functions cj>.j_, the 
potential at a boundary condition control point induced by all 
the singularities on one panel (or subpanel) is 

6 

* = a * 4* 7 + i ^ 1 (ct i' <f> i ) (49) 


The substitution of equation (44) into (49) provides a linear 
relationship between the potential <t> and the unknown doublet 
strength at control points. 

The complete set of linear, algebraic boundary condition 
equations is set up in the following form: 


[A ij ] ^j = RHS i 


(50) 
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where each index i corresponds to the boundary condition equation 

for one control point k' , each index j corresponds to the 

ct'i 

unknown doublet strength at one control point k' , and the 

j 

matrix A-^j and vector RHS^ are arrays of calculated quantities. 
The solution array y j is determined by standard matrix algebra. 


For the perturbation analysis method, it is necessary to 
develop a differential version of equation (50). Recall that the 
panelled configuration geometry is defined by a set of points 
(x,y,z);k (see "Division of a Segment Into Panels"). Now suppose 
that the geometry is distorted by a perturbation to the kth 
point, ( dx, dy , dz ) ^ . The perturbed form of equation (50) is 



The partial derivatives of equation (51) are calculated by follow- 
ing the same procedure used to establish A-^j and RHS^, except 
that each quantity encountered is replaced by its differential. 
For example, consider the differential form of equation (49): 


6 

d<j> = a * d$ _ + <j> -da + E (ct.-d<j>. + <j>.-dct.) (52) 

' ' -i — 1 1 1 1 1 


Each of the differentials on the right-hand side of equation (52) 
has previously been expressed in terms of geometry derivatives 
[equations (45) -(48)] . Consequently, equation (51) reduces to 
the form 


= [B ik 1-dx k + [c ik Wy k + ID ik , ’ dz k 


(53) 


where all the matrix coefficients are calculated quantities. By 
multiplying both sides of equation (53) by the matrix inverse of 
Aji, the desired first-order relationship between geometry pertur- 
bations and singularity distribution is generated: 
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where the partial derivatives in equation (54) are known, calcu- 
lated values. 


Velocity and Pressure Calculation 

The solution array of doublet strengths is converted to sur- 
face velocity and pressure at each panel center, and then the 
forces and moments are integrated. For each body panel, the flow 
velocity in the local panel coordinate system is 


v r 

= V 

* e,. 

+ 

a 


CO 

S 




-> 

->- 



V 

= V 

• e 

+ 

a 

n 

CO 

n 




->• 

-> 



V 

= V 

• e 

+ 

a 

z 

oo 

Z 




(55) 


where a 2 and a 3 are calculated using equation (44). Equation 
(55) expresses the fact that on a body surface the local source 
strength is equal to the perturbation velocity normal component 
and the local doublet strength is the local perturbation poten- 
tial, corresponding to Green's third identity. For steady flow, 
Bernoulli's equation provides the local pressure coefficient. 


S = 1 - 


2 2 ^ 
V r + V + v r 

£ t Z 


V 


ref 


(56) 


Forces and moments are calculated under the approximations 
that the pressure is uniform on each panel and that the flat 
quadrilateral panel is the local geometry. 



GUIDELINES FOR USING THE THREE COMPUTER PROGRAMS 


For the perturbation analysis method, the operation of the 
three computer programs is considerably . more flexible than the 
summary chart of figure 1 suggests. For example, it is not 
necessary to calculate a new partial derivative matrix for each 
angle of attack of interest. Furthermore, the derivative can be 
expressed in terms of standard aerodynamic geometry quantities 
such as thickness and camber distribution, in addition to the 
surface coordinates x, y, and z. The capabilities and operation 
of the three programs is described in greater detail below. 

MCAIR 3-D Subsonic Potential Flow Program (Version 2) 

For a panelled geometry configuration, several flow fields 
can be analyzed simultaneously. Each flow corresponds to a dif- 
ferent specified free stream velocity vector (Voo x , V^y , V ooZ ) and a 

different distribution of prescribed normal velocity (Vjj p ). 
There is one right-hand- side vector in equation (50) for each of 
these flows. 

MCAIR 3-D Geometry Influence Coefficient Program (Version 1) 

Corresponding to any one of the flow field analysis solu- 
tions from the preceding program, this program will generate the 
partial derivative matrix of doublet strength with respect to 
geometry perturbations. Table I presents the geometric quanti- 
ties that are candidates for perturbation. The user identifies 
which of the quantities are to be included in the differentia- 
tion; the others are ignored to save computing effort. 

The structure of the partial derivative matrix is as fol- 
lows. For each control point k^ p , there is one row. Each column 
j corresponds to one independent geometry perturbation, where the 
total number of columns is specified by the user. If dSj is 
defined as the magnitude of the jth differential geometry pertur- 
bation, then the calculated matrix element in column j and row 
k^p is the value of the quantity 

3 u(k<ip ) 

9 Sj 

Each independent geometry perturbation dSj will typically 
represent the change to one of the quantities in Table I. 
However, it is the prerogative of the user to establish a more 
general definition for each dS j , as follows. First, the user 
will designate any number of different quantities in Table I that 
are to be simultaneously dependent on dSj. Second, he will 
prescribe a weighting for each quantity, where that weighting is 
to be the rate of change of the quantity with respect to dSj. 
For example, suppose that the user desires the geometry perturba- 
tion dSj to be the change in thickness at some point on a wing. 
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Table I. Geometric Quantities Subject to Perturbations 


a. Panel Corner Points 
[1 <k<k Q (n s + 1)l 


*1 

x 2 


r 1 

V2 


b. Kutta Condition Direction 
Vectors (One Vector for 
Each Point k£p on a Wake) 


W 


N, 


W N , 


W N , 


x k Yk z k 

111 


Suppose that ZkupR and Z ^-LWR from Table 1 are tbe local upper and 
lower surface coordinates. Then the user would specify a weight- 
ing of +1/2 for dz^ upR and -1/2 for dz-j^^, and the desired 
definition of dSj would be established: 

dz = + 1/2 • dS . 

K UPR 1 


dz = - 1/2 • dS . 

LWR 3 


In addition to the fact that each independent perturbation dSj 
can affect more than one quantity in Table I, it is permissible 
for any single quantity in Table I to be affected by more than 
one independent perturbation dSj. 


dS . 
3 


(dz. 


- dz. 


"UPR 


k LV7R 


MCAIR 3-D Perturbation Analysis Program (Version 1) 

Consider an arbitrary baseline configuration for which a 
panelled geometry representation has been established. Suppose 
that the two preceding computer programs have been executed for 
the panelled geometry at some free stream velocity (V„ x , V^y, V ttZ ) 
and some prescribed normal velocity distribution (Vjsj p ). Then the 

calculations by the two programs can be used to establish an 
input file for this program, the MCAIR 3-D Perturbation Analysis 
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Program (Version 1). For each perturbation case to be analyzed, 
the input to this program consists of an array of geometry per- 
turbation values, dSj. The program will first compute the panel 
corner coordinates of the perturbed geometry and then calculate 
the corresponding surface velocity distribution, surface pressure 
distribution, and forces and moments. 

The only significant difference between this program and the 
first program [MCAIR 3-D Subsonic Potential Flow Program (Version 
2)1 is in the calculation of the doublet strengths In 
this (the third) program, the influence coefficient matrix A-^j is 
not calculated, and the system of linear, algebraic boundary con- 
dition equations is not solved [see equation (50)]. Instead, the 
following first order formula is applied: 


^ k CP } 


*o (k CP> 


I 

j 


( k Qp) 

9 S . 

1 


dS . 
1 


(57) 


where the subscript "o" designates baseline quantities that are 
available in the input file. It is noteworthy that the only 
portion of the perturbation analysis method that is based upon 
first-order approximations with respect to geometry perturba- 
tions is equation (57). Velocity and pressure are not forced to 
vary linearly with respect to dSj, which is expected to be 
significant for wing leading edge radius perturbations. 

Now suppose that one or more additional input files have 
been generated for the baseline configuration, where each file 
corresponds to a different free stream velocity vector an( j a 

different prescribed normal velocity distribution (Vjvjp). 

Identify each file by index IFLOW ( 1 £ I FLOW £ NFLOW ) . In the 
third program, arbitrary linear superposition of flow solutions 
can be constructed from an input set of weightings Wjflow, such 
that the net free stream velocity and prescribed normal velocity 
are 


(Voc , Vco , Vco ) = 

x y z 


NFLOW 
I 

IFLOW 


= 1 [ W 


(Vet 


v a 


Vco ) 


(58) 


IFLOW 



NFLOW 

Z 

IFL0W=1 



IFLOW 


(59) 
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Corresponding to each flow "IFLOW, " equation (57) is used to 
establish y(k£p). Then the net or superimposed solution is 


y(k' p ) = 


NFLOW 
Z 

IFLOW 1 


I, [»-*y] 


IFLOW 


(60) 


The linear superposition associated with the weightings Wiflqw 
does not introduce any error to the perturbation analysis method. 
In other words, the only terms in equations (58) and (59) that 
affect the calculated solution are the net flow conditions on the 
left-hand side. 

An illustration of the utility of equations (58)-(60) is the 
calculation of perturbation analysis solutions at various angles 


of attack 

from 

two input 

files. Suppose 

that 

IFLOW 

= 1 

corresponds 

to 

( Voo X , Vooy , Voo Z ) 

= ( 1 , 0 , 0 ) 

and 

that 

IFLOW 

= 2 

corresponds 

to 

( Voo x * Vooy , Vqo 2 ) 

= (0,0,1). 

In 

each 

case. 

flow 


tangency conditions are assumed (Vjj p = 0). If the angle of 

attack (a) is measured with respect to the x-axis, then the 
weightings for an arbitrary angle of attack are W^ = cos a and W 2 
= sin cl. This procedure was employed in the example solution of 
the next section. 
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EXAMPLE SOLUTIONS 


Two panelled configurations have been analyzed in order to 
evaluate the three computer programs of this report (figure 1). 
The first configuration is a simple wing-fuselage for which the 
sensitivity of a conventional analysis solution to small panel- 
ling gaps was tested. The second configuration is a low aspect 
ratio isolated wing. The perturbation analysis method was used 
to predict the effect of thickness, camber, and twist pertur- 
bations to the isolated wing geometry. All the calculations of 
this section were performed on the McDonnell Douglas CDC CYBER 
750 computer. 

Conventional Analysis of a Wing-Fuselage 

Figure 18 is a sketch of two panelled geometry represen- 
tations of a fighter-type wing mounted on an axisymmetric fuse- 
lage. In the panelling of figure 18a, there are no unpanelled 
gaps on the configuration surface. The panelling of figure 18b 
is identical, except that approximately one-half of the stream- 
wise defining points on the fuselage were removed and the 
vertical plane of the wing tip was not panelled. The reduction 
in fuselage panel density introduced small vertical gaps between 
the fuselage and wing panel edges at the exposed root, where the 
size of the gaps is on the order of one-half of one percent of 

the local wing chord. The objective is to determine whether the 
MCAIR 3-D Subsonic Potential Flow Program (Version 2) can gener- 
ate an accurate pressure distribution for the geometry with gaps 
in the panelling (figure 18b). 



a) No Gaps in Paneling b) Gaps at Wing Root and Tip 

Figure 18. Two Paneled Representations of a Wing-Fuselage 
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As a standard of comparison, the PAN AIR pilot computer pro- 
qram of reference 3 was applied to the densely paneled geometry 
of figure 18a. The two calculated pressure distributions are 
nearly identical, as illustrated in the vicinity of the wing root 
by figure 19. It is concluded that elimination of all gaps in 
the panelled geometry is not necessary. The calculations for 
this example were performed as part of the 1980 MCAIR Independent 
Research and Development Program. 

Perturbation Analysis Method for an Isolated Wing 

The wing panelling of figure 20 was selected as the baseline 
geometry for demonstrating the perturbation analysis method. At 
each span station, the wing section is an untwisted NACA 0012 air- 
foil (12% thickness ratio). Each of the five perturbed wing 
geometries of figure 21 was analyzed by both the conventional 
panel method [MCAIR 3-D Subsonic Potential Flow Program (Version 
2)] and the perturbation analysis method. A comparison of the 
calculations will establish an indication of the accuracy of the 
perturbation analysis method. 


Each of the perturbed wing geometries was generated by coor- 
dinate displacements normal to the baseline wing chord plane 
(Az) . No perturbations in the x or y-directions were considered. 
Ruled surfaces were used to define the panelling between the root 
and tip sections, with one exception. For Wing B, the twist dis- 
tribution was specified to vary linearly across the span. All of 
the wings (except E) were analyzed at two angles of attack, 0° 
and 5°. Here, angle of attack is measured with respect to the 
baseline wing chord plane. 

The conventional panel method solution for the baseline wing 
is presented in figures 22-25. For verification, a solution by 
the PAN AIR pilot computer program is also presented. For PAN 
AIR only, additional panels were distributed across the open wing 
tip of the baseline geometry to eliminate unpanelled gaps. The 
section properties c £, C(j, c m of figures 23 and 25 are based upon 
the local chord, and the moment reference point for c m is the 
quarter chord of the local wing section. It is noteworthy that 
at zero incidence the inboard and outboard section drag 
coefficients are approximately +100 and -100 counts, respectively 
(figure 23). However, the total wing drag coefficient is less 
than one count for both panel methods, which is consistent with 
the potential flow law that nonlifting bodies have zero drag. 

Perturbation analysis solutions for wings A through E at 0° 
and 5° angle of attack are presented in figures 26-43. The 
legend "exact solution" in the figures designates the 
conventional panel method solution for the perturbed wing 
geometry. It is apparent that the perturbation analysis method 
is quite accurate for a variety of large geometry perturbations. 
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Wing Root Pressure Distribution 
(35% Semispan) 



Figure 19. Calculated Pressures Near Wing Root 

4° Angle-of-Attack 







Figure 20. Baseline Wing Paneling for Demonstration of Perturbation Analysis Method 
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Baseline Wing 
(NACA 0012) 


Wing A 
(NACA 0006) 


Wing B 
(10° Twist) 



Perturbation Analysis Method 











Unlike the existing small disturbance (linearized) methods, the 
accuracy does not diminish substantially in the leading edge 
region, even if the radius of curvature is changed. 

The only appreciable difference between the perturbation and 
exact solutions is the pressure distribution on the lower surface 
of the wing with deflected flaps and the aft lower surface of the 
supercritical wing. The explanation for this difference is 
apparent upon examining figures 22 and 42. The 20° flap deflec- 
tion changes the level of lower surface pressure coefficient from 
the baseline level Cp a, 0 to Cp r ul/2. This represents one-half of 
the theoretical maximum perturbation (to stagnation). The pertur- 
bation analysis method does not account for the fact that the 
lower surface velocity cannot decrease indefinitely with respect 
to increasing flap deflection; therefore, the perturbation analy- 
sis method slightly overestimates the magnitude of the local 
pressure coefficient. 

Based on figures 26-43, it is concluded that for section 
geometry perturbations representative of preliminary aerodynamic 
design of clean wings, the accuracy of the perturbation analysis 
method is competitive with conventional panel methods. The 
advantage of the perturbation analysis method is the small 
computing effort required for each additional perturbation case. 
In the present example, each perturbation case required only 7 
seconds computing time versus 245 seconds for the conventional 
analysis solution. The generation of the partial derivative 
matrix by the MCAIR 3— D Geometry Influence Coefficient Program 
(Version 1) required 925 seconds each for 0° and 90° angle of 
attack. In the present example, it is cost efficient to employ 
the perturbation analysis method if nine or more cases are to be 
analyzed. 
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Figure 22. Baseline Pressure Distribution at 
0° Angle-of-Attack 

NACA 0012 
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Inboard 
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Figure 28. Pressure Distribution of Wing A at 5° Angle-of-Attack 

NACA 0006 
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Figure 30. Pressure Distribution of Wing B at 0° Angle-of-Attack 

10° Twist 
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Outboard 
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Figure 40. Pressure Distribution of Wing D 
at 5° Angle-of-Attack 




Exact solution 0.318 

Perturbation 
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Figure 42. Pressure Distribution of Wing E at 0° Angle-of-Attack 

20° Flaps 




0 



0 



Figure 43. Force and Moment Distribution of Wing E at 0° Angleof-Attack 

20° Flaps 
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CONCLUSIONS 


The perturbation analysis method is similar to classical 
small disturbance methods in the sense that linearized boundary 
conditions are employed. However, the perturbation analysis 
method is significantly more accurate for two reasons. First, 
the linearization is performed on an arbitrary baseline geometry, 
not a flat plate. Second, only the doublet strengths at discrete 
points are linearized, not surface velocity and pressure. 

A.fter the partial derivative matrix has been established for 
a baseline configuration, each perturbation analysis is more than 
an order of magnitude more efficient than a conventional panel 
method solution. The reason is that there is no aerodynamic 
influence coefficient matrix to calculate nor is there a large 
system of linear, algebraic boundary condition equations to 
solve. It has been demonstrated that the accuracy of the pertur- 
bation analysis method is competitive with conventional surface 
panel methods for analyzing the effect of large changes in wing 
section geometry. 

The utility of the perturbation analysis method for prelimi- 
nary design is obvious. There are two additional applications 
for which the method is expected to be equally useful. The first 
is in iterative solutions to prescribed pressure design problems, 
where the cost of analyzing the updated geometry of each itera- 
tion cycle would no longer be a deterrent. The second applica- 
tion is the prediction of strong viscous-inviscid interactions by 
iteration between viscous and inviscid flow analyses. The number 
of iteration cycles and/or angles of attack would not be limited 
by the high cost of conventional inviscid flow analysis. 
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APPENDIX A 


QUADRATIC SURFACE FIT THROUGH 
TWO SETS OF THREE COLINEAR POINTS 


Consider the points & (1 <_l<_ MAX) in figure A-l, where MAX 
is either 6, 7, 8 or 9. The set of points [1,2,4] is colinear, 

and the set of points [1,3,5] is colinear. The objective is to 
establish a quadratic function y(5,n) that will match arbitrary 
values y^ exactly at points 1 <^£_< 5 and approximately at points 
6 £&j£MAX. The given (known) quantities are 0, s jj, for 2 <_Z<_ 5, 
MAX, and ( n) g for 6 MAX. Any or all of the values - S£ 

(2 <_£< 5) can be negative or positive. Assume that the function 
y(£,h) is expressed in the form 


y(?,n) = u 1 + + 3 2 r ' + 8 3 S 2 + 3 4 Cn + B 5 (n 2 -? 2 tan 2 e ) (Al) 


v 



Figure A-1. Nomenclature for Two Sets of Three Colinear Points 
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Then the coefficients 3i can be readily established from the 
following formulas: 


6^* 2cos0 


s 5 2(u 3" u 1 ) " s 3 2 (u 5 -p 1 ) 

, X < 

s 4 2 ( p 2 - Ui ) - s 2 2 (u 4 - Pl ) 

8 3 8 5 (s 5- s 3 ) 

T * 

s 2 s 4 ( s 4 -s 2 ) 





(A2) 


B 2 • 2 sin 0 


[ S 5 2 (y 3 ' y l ) ' s 3 2 (w 5 _,J 1 ) ' 

[ S 4 2 (p 2 ~ v 1 ) ~ S 2 2 (u 4 ~ u i y 

1 S 3 S 5 (S 5 - S 3 ) 

t S 2 S 4 (S 4 - S 2 > 


(A3) 


B 3 *2cos 6 


S 3 (y 5 _y l > " S 5 (y 3 _y l > 

' S 2 (y 4 _y l ) " s 4 (y 2 " y l ) ' 

S 3 S 5 (S 5 " S 3 ) J 

h s 2 s 4 (s 4 -s 2 ) 


( A4 ) 


B . " 2 cose sine 

4 


' s 3 (y 5 ' y l ) " S 5 (y 3 " y l ) ] 


s 2 ( y 4 _y l) s 4^ y 2 y l^ 

S 3 S 5 (S 5 " S 3 ) j 


S 2 S 4 (S 4 " S 2 ) 


(A5) 


MAX 

l SGN , 
2 = 6 


{ y * ~ “j. - V 


~ B 3 5 2 " 6 4 5 £. n £j 


MAX 299 
£ |n t - e • tan^e | 
2 = 6 *■ 


( A6 ) 


where 


SGN 2 


-1 if 2 = 6 or 8 
+1 if 2 = 7 or 9 
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The above definition of terra SGN^ is based on the assumption that 
the quantity (r | 2 - ?£ tan 2 0 ) is j 0 if £ is 6 or 8 and > 0 if £ 
is 7 or 9. This assumption is always valid for the panel method 
formulation of this report. 

The orientation of the (E , , n ) -coordinate system of figure A-l 
and the functional form of equation (Al) were selected to 
simplify the determination of equations (A2) through (A 6 ). The 
derivation is straightforward but will not be presented here 
because a simple substitution into equation (Al) validates the 
formulas for 3j_. 

It is noteworthy that each 3^ is linear with respect to y^ . 
Specifically, equations (A2) through (A 6 ) are in the form 




> 


6 . c 
i5 


MAX 
E (b 
i=6 




(1 £ i 1 5) (a7 ) 


where 6-^5 = (0,1) if i is (^5, =5), and where the coefficients bj.£ 
are explicit functions of the given quantities defined earlier. 

Now suppose that the given quantities include known values 
for (1 _< £_< MAX) and that it is desired to determine 33i/30, 

3 3j_/ 3s £ for 2 <£< 5, and ( 3 B±/ 3€£, 3 Bj./ 3n£) for 6 ££< MAX. This 
is accomplished by differentiating equations (A2) through (A 6 ). 
First, it is convenient to introduce twelve quantities defined by 
equations (A 8 ) and (A9). 


\ “ S 5 (y 3 y l ) 

- s 3 2 (p 5 -m 1 ) 

A D = S 3 S 5 (S 5‘ S 3 ) 


B N E S 4 2(y 2" y l ) 

- s 2 2 (y 4 -p 1 ) 

B D E S 2 S 4 ^ S 4 _S 2 ^ 

(AS) 

C N E S 3 (y 5 _y l ) 

- s 5 (y 3 -y 1 ) 

C D E S 3 S 5 (S 5 -S 3 } 


D N E S 2 (y 4“ y l ) 

- s 4 (y 2 - P;L ) 

D D E S 2 S 4 (S 4“ S 2 ) 
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f a = <VV 

F c 5 <V C D> 


f b 5 <VV 


F D 5 <VV 


(A9) 


The partial derivatives of F A , F B , F c , and F D follow 


3f a _ 

3 S 2 

{ 'V S 3 (y 5" V l ) 

' A N S 5 (S 5 2S 3 ) }/A D 


8F a _ 

3 s c 
5 

{+A D 2s 5 (y 3" y l } 

+ a n s 3 (s 3~ 2s 5 )}/A D 2 


3F b _ 
3s 2 

{ 'V S 2 (u 4" y l ) 

' B n S 4 (S 4 -2S 2 ) }/B D 

(A10) 

9F b _ 

3S. 

4 

{+B D 2s 4 ^ y 2 _y 1^ 

+ B N S 2 {S 2" 2S 4 )>/B D 2 


3F c _ 

3S 3 

{+ c D (y 5 - y l ) ' 

C N S 5 (s 5 _2S 3 )}/C D 2 


8F C 

3S 5 

{ - c D (y 3 - y i ) + 

C N S 3 (S 3" 2s 5 ) }/C D 


8f d _ 
9s 2 

{+D D (y 4” y l ) ' 

°N S 4 ( s 4 -2s 2^ }/D D 


3F d _ 

3s 4 

{ ' D D (y 2" y l ) + 

°N S 2 (s 2~ 2s 4 ) }/D D 
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With the aid of the preceding quantities, equations (A2) 
through (A6) are differentiated to establish equations (All) 
through (A16). 


2cos0dg i 


2g^sin0de + 


+ 



(All) 


2sin0d02 = -202cos0d8 - 


2cos SdB^ = 4 B 2 c os 0 sine d 0 


+ 


+ 


2cos0sin9dB 4 = -2B ^cos ( 20 ) d 0 - 



(A12) 


(A13) 


( A14 ) 
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dg 5 ( TDEN^ 


MAX 
I SGN 
£ = 6 


£ 


S £ dB l n £ dB 2 ” ? £ dB 3 " 5 £ n £ d3 4 
+ ( " B 1 - 2? £ B 3 - V4 )d? £ 


+ ( “ 6 2 - ^4 )dy £ 


+ ( 


TNUM 


MAX 

* 

2 1 

IK 2 sine^ 

\ I 

Z ( 2 • S GN ) • . 

Q — (Z 

S £ tan ed? £ - n £ dn £ + 

I £ 

de 

3 Q I 

X/ — o 

- 

l cos 0 } 

J 


where 


(A15 ) 


MAX 

TNUM 5 Z SGN 
£ = 6 


y ° ~ y-1 B l^£ “ B 2 y £ 


6 3^£ “ B 4 ? £ y £ 


( A16 ) 


MAX „ 

TDEN = z |n/ - 5 • tan e 

£ = 6 X 


The desired partial derivatives of 3 i are expressed 
implicitly by the coefficients of the differentials in equations 
(All) through (A15). 
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APPENDIX B 


DIVISION OF A PANEL INTO SUBPANELS 


The division of an arbitrary noncoplanar quadrilateral panel 
into eight triangular subpanels is based upon the approach of 
reference 3. For each subpanel, a quadratic doublet distribu- 
tion is established by a fit to the doublet strength at several 
spline control points on the panel. A sensitivity matrix is 
generated for which each matrix element is the contribution of a 
unit strength doublet density at one spline control point to one 
of the six coefficients of the quadratic distribution. The 

present approach differs from that of reference 3 in the sense 
that the sensitivity matrix is generated from simple, explicit 
mathematical functions and no matrix inversion is required. The 
mathematical formulation is presented in this Appendix. Also 
presented is the procedure for generating a second sensitivity 
matrix, one that relates the quadratic doublet coefficients to 
panel geometry perturbations. 

The division of a noncoplaner panel into eight triangular 
subpanels is illustrated in figure B-l. The subpanel vertices of 
figure B-l are the same nine points 1 < k < 9 that were depicted in 
figure 9. Similarly, the coordinate displacement arrays X-j^ 
defined by equations (3) and (4) also apply to this Appendix. It 
is easily proved that the four interior subpanels (5 < i sp <8) are 
coplanar, regardless of whether the four corner points (k = 
1,2, 3,4) are coplanar. 
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The nine points k are the spline control points. Suppose 

that an arbitrary value of doublet strength (y-^) is assigned to 
each spline control point. For each subpanel, it is desired to 
establish the quadratic doublet distribution that will satisfy 
the following conditions: 
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(i) Consider the three spline control points at the sub- 
panel vertices and the two spline control points that 
lie on extensions of the subpanel edges. At each of 
these five points k, the quadratic distribution is to 
match . 

(ii) Consider the subpanel edge that separates interior from 
exterior subpanels. At the midpoint of this edge, both 
the doublet strength and the derivative of the doublet 
distribution in the direction normal to the edge are to 
match the corresponding values for the adjacent 
subpanel. This will generate continuity and slope 
continuity at the edge midpoint. 

The above two conditions imply that for each pair of adja- 
cent subpanels, the quadratic doublet distributions will match at 
three points along the common edge (or an extension of that 
edge) . This matching is (usually) enforced even if the adjacent 
subpanels are not part of the same full panel. Along an edge, 
the quadratic surface distribution is a parabola. Since there is 
only one parabola that can pass through a set of three points, 
the doublet distribution will be continuous across the entire 
common edge of each pair of adjacent subpanels. 

Figure B-2 illustrates a local coordinate system ( £, n, C) 
that is defined for each subpanel igp. The C -axis is perpendicu- 
lar to the subpanel, and the origin is at the one subpanel vertex 
that is not coincident with an edge midpoint of the full panel. 
The index & denotes the spline control points in an order differ- 
ent from index k (see Table B-I ) . The ordering of £ was selected 
to be consistent with the method of Appendix A, which will be 
employed in establishing the quadratic doublet distribution. For 
each subpanel igp, "the nomenclature k(£) and £(k) refers to the 
conversion between £ and k that is provided by Table B-I. 


n 
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Figure B-2. Local Coordinate System for a Subpanel 



Table B-l. Conversion Between Index k and Indices ( I , igp) 


k (for 1 < i sp < 4) 

k (for 5 < igp < 8) 

'SP 

9 

(i S p + 4) 

(igp - 1) if igp #5 
8 if igp = 5 

(igp + 3) if igp ¥= 1 
8 if igp = 1 

'SP 

(igp + 1) if igp =£ 4 
1 if igp = 4 

(igp + 1 ) if igp ^ 8 
5 if igp = 8 

(igp - 1) if igp =£ 1 
4 if igp = 1 

(igp + 2) if igp < 6 
(igp — 2) if igp > 6 

9 

(igp - 4) 

(igp + 5) if igp =£ 4 
5 if igp = 4 

(igp — 3) if igp 8 
1 if igp = 8 

(igp + 6) if igp < 2 
(igp + 2) if igp > 2 

(igp — 5) if igp # 5 
4 if igp = 5 

(igp + 2) if igp < 2 
(igp - 2) if igp > 2 

(igp — 2) if igp ^ 6 
(igp —6) if igp > 6 


Quantities designated {P ± £, 9, s £ , j , E, 2 > r] 2> h • n 3 » n M i are 
to be calculated for each subpanel i g p. The formulas are 
provided below: 


p = p = y 

u ' a " i,kU) 


X i ,k (1) 


(1,1) £ (i,£) £ (3,5) (Bl) 


0 = k Cos -1 [ (P -P „)/( IpJ • IpJ ) ] 
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^2 = s 2 cos0 
n 2 = _s 2 sine 

£ 3 = s 3 cose 

n 3 = sine 


(B5) 


? M 2 ^2 + ^3* 


(B6) 


n M = 2 (n 2 + n 3 } 


J 


The method of Appendix A is to be employed in order to estab- 
lish the quadratic function that will exactly match yi c (£) for 
1 < £ < 5 and also match some dummy value y^ at the point 
where subscript "M" designates the midpoint between sub- 
panel vertices 2 and 3. The required given values for Appendix A 
are 9 from equation (B2), S£ from equation (B3), “MAX" which is 
6, and the coordinates ( ?M' ’"Im) from equation (B6). It is note- 
worthy that the coordinates designated here correspond to 

(56/^6) in Appendix A. From these given values, the method of 
Appendix A will provide an array of coefficients bj_£ 

[(1,1) < (i,£) (5,5)] and a single coefficient b 5 M such that 


yU/n) = y k(1) + 3 ] _C + S 2 n + 8 3 ? 2 + 8 4 ?n + 8 5 (n 2 ~Z 2 tan 2 e ) (b7) 


where 


8 . 


l 


5 

£ = 1 b i^ ,W k(£) + 6 i5 b 5M* U M 


( B8 ) 
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Regardless of the value for equations (B7) and (B 8 ) 

satisfy condition (i) . It remains to determine y M such that 
condition (ii) is satisfied. For doublet continuity, the value 
of y^ f° r eac h exterior subpanel igp must be equal to the value 
of y m for the adjacent interior subpanel (igp+4). For slope 
continuity, 

{vw(« M ,n M )-I+(n 3 -n 2 )S 5 - < 6 3 “S 2 } ] }-i sp 
= {»vU M .l M > • [-(n 3 -n 2 )e 5 + <5 3 ' e 2 ) e n '} (i gp +4) (B9) 


With the aid of equations (B7) and (B 8 ), equation (B9) can be 
expressed in the following form: 


Pm = ^l{ Cr[VkU)li SP + d *' ty *(M ] (i S p+ 4 )} 


where C£ and d£ are calculated coefficients that are independent 
of y k . 


By substituting equation (BIO) into (B8) and then redefining 
the coefficients b^£ to include the contributions of C£ and d£, 
equation (B 8 ) becomes 


B i = 


5 

E 

£=1 


{b i£ 


' W k(£) } 


6 .- 

i5 


8 

E {b 
£ = 6 


5£ * y k ( ft ) 


> 


(Bll) 


For each subpanel igp (l£igp£8), equations (B7) and (Bll) 
provide the desired quadratic doublet distribution. Each coeffi- 
cient (l£i£5) is a linear function of y k (l£k_<9). The 

sensitivity matrix elements bj_£ [(1,1) £ (i,ft) £ (5,5)] and b 5 £ 
(6 £ft£ 8 ) are calculated independently of y k by applying the 
formulation described above to the nine given quantities Xp k , 
X 2 k , X 3 k (2£k£4); recall that Xu = X 21 = X 31 = 0. If one of 
the edges of the full panel is zero length, the above formulation 
requires only two slight modifications. First, the two exterior 
subpanels that lie on the zero length edge collapse to zero area 
and are ignored. Second, for each of the two interior subpanels 
that touch the zero length panel edge, the linear relationship 
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between and yj. is established by passing a parabola through y^ 
at the three spline control points on the panel edge that touches 
the subpanel point 

The preceding portion of this Appendix described the proce- 
dure for establishing the quadratic doublet distribution corre- 
sponding to known values X-j/k and arbitrary unknown values y^. 
Now suppose that each y k is assigned a fixed, known value. The 
method for determining the first-order change in quadratic 
doublet distribution induced by geometry perturbations (dX^) is 
presented below. 

Define the array v m (l<ra£9) by equation (18). Then dXi k 
for l£k£9 can be expressed by equations (19) and (20). For 
each subpanel igp, the derivative of equation (Bl) provides the 
formula for dP^j,: 


is. “ jf x { DPI Hj' dv [i+3(j-l) ]} (1 ' 1 ^ 1 < i .*■) 1 (3,5) (B12) 


where DPI £j = {DXI k(l)f . - DXI k(1)/j > 


(B13) 


By applying the calculated values DXI k j and DPI to the deriva- 
tives of equations (B2)-(B6), the following quantities are calcu- 
lated for each subpanel igp: 


30 

3S * 

3a ij 

3e 2 

3n 2 

3^ 

3n 3 

3£ 

M 

M 

3 V ' 
m 

3 v ' 
m 

3v ' 
m 

3 v ' 
m 

3 V ' 

m 

3 v ' 
m 

3 V ' 
m 

3 v ' 
m 

3 v 

m 


(1,1, 1,1) < (i,j,£,m) < (3, 3,5,9) 


87 



Calculated values for 33-j_/36, 33i/3s^, 33^/3^ M , and 33-^/3n M 
are generated by applying the procedure of Appendix A to the fol- 
lowing known values for each subpanel igp: 

{ 8 / y k { £) (1 5 )- 


The value y m = is obtained from equations (B7) and 

( Bll ) • Again, "MAX" is 6 and the index &=6 of Appendix A corre- 
sponds to the midpoint "M" here. The formulation of Appendix A 
is based on the presumption that is a fixed quantity 

(dP M = 0). Here, however, is dependent upon geometry pertur- 

bations (dy M ^ 0), and the differential expression for d3^ is 
written in the form 


d3 ± 


33 . 

l 

39 


dO 



33. 


3C 


az 


M 


M 


33 ± 

+ — — dri 


3n 


M 


R + 6 i5 b 5M* dlJ M 


( B14 ) 


All the coefficients of equation (B14) are known quantities, 
including bg^ which appeared earlier in equation (B8). It 
remains to establish formulas for 33i/3v m . 



+ 


6 i5‘ b 5M' dy M 


(B15) 


9 

E 

m=l 


(DB . • dv ) 

im m 


S i5 ' b 5M‘ dy M 


In order to determine the relationship between dyj^ and dv m , 
equation (B9) is written in the following form: 


* 

5 


► * 
5 

■ ^ (g.-e.) 

i=l 1 1 

- + . 
1 SP 

Z (g. 8.) - 

i=l 1 1 

» < 


(i SP +4) 


= 0 


( B16 ) 
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where 


g l = 71 3" 11 2 
g 2 = C 2~ ? 3 
g 3 = 2 ? m’ 3 -n 2^ 

g 4 = [n M (n 3 -n 2 ) + £ M (S 2 ~ S 3 )] 

2 

g 5 = 2[n M (5 2 -5 3 ) " tan 


The substitution of equation (B15) into the differential form of 
equation (B16) results in the following equation: 


9 

' (g 5' b 5M )dl, M + S , 

[\\ 

i 3g i\ 

i g i ' DB im + ®i av j 

dv • 
m 

m=l 


\ m/ J 

* 


( B18 ) 



9 

r 5 / 

3g.- \ ’ 

■ 

+ ■ 

(g 5' b 5M )dM M + = 

Z I g. ■ DB . + e . 

. lm l 3v / 

dv 

m 


L m=l 

Ll=l\ 

m/ j 



(i SP +4) 


0 


Equation (B18) can be written as 


9 

dy = Z (DC -dv ) 
M _ , mm 
m=l 


( B19 ) 


where DC m is calculated by consolidating the coefficients of 
equation (B18). Finally, when equation (B19) is substituted into 
equation (B15), the accumulated coefficients of dv m are the 
desired values 38i/3v m [ ( 1, 1 ) £ (i,m) £ (5 , 9 ) 3 . If one of the 

edges of the full panel is zero length, the procedure for 
calculating 3 8^/ 3v m requires two modifications that correspond to 
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those described earlier. First, the two exterior subpanels that 
lie on the zero length edge are ignored. Second, for each of the 
two interior subpanels that touch the zero length panel edge, 

dy M = °* 


Now consider arbitrary simultaneous perturbations to both 
panel geometry (dv m ) and doublet strength at the spline control 
points (dy k ); see equation (Bll). For each subpanel igp, the 
change in quadratic doublet coefficient (d3i> can be predicted by 
the following formula: 


dB 


9 

E 

m=l 



+ 


5 

Si l 1 b i!i' dy kU) 


+ <5 . _ 
i5 




(B20) 


All the coefficients of equation (B20) are known quantities, 
calculated by the procedures presented above in this Appendix. 
In the same sense that the linear relationship between 3-j_ and 
of equation (Bll) is appropriate for a conventional panel method, 
the linear relationship between d3i and dv m of equation (B20) is 
appropriate for the analysis of small geometry perturbations. 
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